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Abstract 

The  area  between  the  continental  shelf  and  the  Gulf  Stream  is  known  as  the  slopewater 
region.  In  the  past  fifty  years  several  experiments  and  studies  of  this  region  have  taken 
place  with  the  Mid- Atlantic  Continental  Slope  and  Rise  (MASAR)  experiment  being  one 
of  the  most  recent.  Csanady  and  Hamilton  (1988)  compiled  all  the  known  information 
and  data  from  the  slopewater  region  and  developed  a  simple  dynamical  model  of  the 
flow.  Based  on  this  model’s  transport  stream  function  and  Stommel’s  Gulf  Stream 
model,  finite  centered  differencing  was  used  to  develop  a  numerical  scheme  of  slopewater 
circulation. 

The  model  was  first  developed  using  Stommel’s  parameters  for  circulation  within 
the  North  Atlantic  Gyre.  Stommel’s  model  was  used  as  the  basis  for  the  new  scheme 
in  order  to  calibrate  the  model  with  his  exact  solution  of  the  streamfunction  for  the 
North  Atlantic  Gyre.  Once  verified,  Stommel’s  parameters  were  replaced  by  Csanady 
and  Hamilton’s  values  for  slopewater.  This  is  a  report  on  the  development  of  a  new 
numerical  model.  It  is  also  a  comparison  of  the  new  scheme  to  both  Csanady  and 
Hamilton’s  model  and  an  observational  schematic  for  the  region  from  the  MASAR 
experiment. 

Keywords:  Circulation-Slopewater-Modeling-Physical  Oceanography 
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Circulation-  A  Short  History 


Ocean  circulation  can  be  described  as  either  wind-driven  or  thermohalinp  Th 

foil  wing  conditions:  m  1947.  He  assumes  the 

•  the  pressure  gradient  is  balanced  solely  by  the  Coriolis  force, 

•  horizontal  velocities  and  pressure  gradient  disappear  with  depth, 

•  lateral  stresses  are  neglected, 

•  and  the  flow  is  steady-state. 

The  Westward  Mensifi. 

erations  that  could  develoo  within  th  1  r  er  to  counterbalance  any  accel- 

forces  in  the  ^^0'  wTnd  st^s  and  b  n  “  ?■  ^  m°deI’  St°mmel  added  Actional 
to  develop  a  general  solution  to  ocean  °  ^  ”Ctl°n'  ^ese  assumptions  allowed  him 
1950,  Munk  w rZ  on  tZZ  TZtT"  T*""  Tw°  *eare  la‘« 

bines  the  work  of  Sverdup  and  S  tom  "m  T  I"  hiS  reSCarch-  M"ak 

ocean  circulation.  He  a  baZ  ”  /I  *  r^^taiioa  of 

at  a  slight  angle  to  each  other  rather 

force,  a  basin  that  encompasses  both  nrirf.u  ^  ,  '  , eraJ  lnct,on  as  the  dissipative 

observed  wind  instead  of  a  sinusoidal  function.  ^  S°Ut  ^  hemispheres’  and  uses  the 

Small  contributions  were  made  to  the  understanding  of  the  Gulf  St™™  d  , 
water  region  over  the  next  thirty  years.  The  next  lar  /  l,  Stream  and  sloPe- 
recently.  In  the  mid-80’s  the  MiH  Ati  *•  e  •  ge  contribution  was  added  more 

experiment  washed  ^  TWs'  p^m  T*  °C“"  Hl“  (MASA*> 

Mineral  Management  Service  fa  U  S  Iv  WaS  designed  to  meet  the  objectives  of  the 
6  Service  (a  U-S-  government  department)  which  were 

•  to  determine  the  general  circulation  features  on  the  continental  slope  and  rise 
’  ir  “ty  °f  a^  other 

— — — — * 


The  MASAR  experiment  involved  several  cruises  in  the  region  to  put  into  place  many 
data  collection  tools  such  as  current  meters  and  to  collect  water  samples  and  real  time 
data.  Csanady  was  the  key  person  for  the  interpretation  of  data  for  upper  slope  circu¬ 
lation  in  the  MASAR  experiment.  He  compiled  all  the  major  sources  of  information  on 
the  slopewater  region  and  produced  a  simple,  linear  dynamical  model  of  circulation  in 
the  area.  His  model  is  developed  under  the  following  assumptions: 

•  the  flow  is  steady-state, 

•  Coriolis  is  balanced  by  wind  stress, 

•  and  the  region  is  rectangular. 

The  model  also  includes  non-zero  boundary  conditions  for  the  streamfunction.  Csanady 
and  Hamilton’s  paper  Circulation  of  Slopewater  (1988)  included  an  observational  sce- 
matic  (Figure  1)  from  the  MASAR  experiment  showing  the  details  of  Slopewater  circu¬ 
lation  as  it  thought  to  occur.  This  paper  is  the  springboard  for  the  numerical  scheme 
that  was  developed  to  reproduce  this  schematic. 


2  The  Model 

The  development  of  the  model  began  with  the  following  basic  assumptions: 

•  the  water  parcel  is  a  homogeneous  layer  of  incompressible  fluid, 

•  the  region  is  enclosed  in  a  rectangular  box, 

•  there  are  zero  boundary  conditions, 

•  and  the  flow  is  steady-state. 

The  density  is  assumed  to  be  constant  throughout  the  parcel  because  it  is  a  relatively 
shallow  layer  of  approximately  500  meters.  Therefore  it  is  reasonable  to  assume  incom¬ 
pressibility.  The  formulations  that  are  the  foundation  of  modeling  slopewater  circulation 
are  presented  in  sections  2.1  through  2.4  of  this  paper. 

2.1  Compressibility 
The  compressibility  of  a  liquid  /?  is  defined  by 

0  =  -L*!L  [lmdV\{dt\  fldV\/dp\ 

V  dp  \V  dt)[dp)  [ V  dt)  •  \dt)  W 

where  /?  is  the  compressibility,  V  is  the  volume,  p  is  the  pressure,  %  is  the  change  of 
volume  with  time,  and  ^  is  the  change  of  pressure  with  time. 

If  the  fluid  is  incompressible  then  (3  must  be  zero.  Therefore 
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If  the  mass  of  fluid  is  constant  over  volume,  V  and  density  is  p= m/V  then 


1  ff£  _  Y_  £  (™\  _  _1  <w_ 

p  dt  m  dt  \V )  V  dt 


(3) 


Hence  incompressibility  occurs  when  f’^"  =  0orp^t=0  [Pond  and  Pickard,  1983]  . 
Incompressibility  is  an  important  assumption  in  this  problem  for  it  allows  the  stream 
function,  ip,  to  exist. 


2.2  Stream  Function 

The  total  material  derivative  of  the  fluid  density  is  known  to  be 

Dp  dp  _  _ 

Dt  dtT  p 

For  an  incompressible  fluid,  the  density  does  not  change.  In  this  case 

=  o  • 

Dt 


The  equation  of  continuity  which  is  given  by 

Tt +  ^  = 0 

can  be  expanded  in  the  form  of 

dp  „  _  _  _ 

at 

By  substitution,  the  equation  of  continuity  of  an  incompressible  fluid  becomes 


(4) 

(5) 

(6) 

(7) 


V  •  v  =  0, 


(8) 


which  is  expressed  below  for  a  two-dimensional  flow  velocity  ,  v  =  (uyv)  as 


du  dv 
dx  dy 


=  0 


(9) 


For  the  equation  of  continuity  to  be  satisfied,  let  u  =  and  v  =  —  where  rj)  is  the 
two-dimensional  streamfunction  of  fluid  dynamics.  The  velocity  field  v  is  now  expressed 
in  vector  form 

.  M 

When  ip  is  a  constant,  the  streamfunction  produces  curves  known  as  streamlines.  A 
streamline  is  a  curve  in  space  drawn  so  that  the  velocity  vectors  are  tangent  to  the 
curve.  This  allows  the  velocity  field  to  be  represented  by  streamlines. 
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2.3  Stommel’s  Gulf  Stream 

Stommel’s  Gulf  Stream  model  is  the  pattern  for  the  slopewater  model.  His  model  is 
based  on  the  following  assumptions: 

•  the  ocean  is  rectangular  with  the  y-axis  pointing  northward  and  the  x-axis  east¬ 
ward 

•  the  boundaries  of  the  rectangle  are  between  0  <  x  <  X  and  0  <  y  <  b 

•  the  ocean  is  a  homogeneous  layer  of  constant  depth  D  when  at  rest 

•  when  there  are  currents,  the  total  depth  is  D  +  h  where  the  depth  varies  by  h 
which  is  much  smaller  than  D 

•  functional  form  of  the  wind  stress  acting  over  the  area  is  -F  cos(^) 

•  the  frictional  dissipative  forces  are  -Ru  and  -Rv  where  R  is  the  coefficient  of 
friction  and  u  and  v  are  the  velocity  components 

•  the  Coriolis  parameter  f  is  introduced  as  a  function  of  y. 

The  steady-state  equations  of  motion  are  written  in  the  following  manner: 

The  equation  for  motion  in  the  y  direction  is 

0  =  f(D  +  h)v  -  F cos( ““ )  -  Ru  -  g(D  +  h)j£  (11) 

The  equation  for  motion  in  the  x  direction  is 


0  =  -f(D  +  h)u-Rv-g(D  +  h)j±.  (12) 

These  equations  are  cross-differentiated  and  the  equation  of  continuity  is  applied  to  the 
equations  of  motion  resulting  in 

To  a  first  approximation,  h  is  negligible  compared  to  D.  This  allows  the  previous  equa¬ 
tion  to  be  rewritten  as  4 


ctv  +  7  sin 


dv  du 

+  75 - 7—  =  0 

ox  dy 


(14) 


where  a  =  (£)  (§£)  and  7  =  &. 

The  ocean  is  considered  to  be  incompressible  which  allows  the  introduction  of  the  stream 
function.  From  the  previous  equation,  the  general  equation  is  stated  as 


(15) 
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and  assumed  boundary  conditions  are 

=  0(A,y)  =  ip(x,0)  =  ip(x,b)  =  0  .  (16) 

Stommel  solves  the  nonhomogeneous  differential  equation  by  inspection  and  separation 
of  variables  for  ip  (See  Appendix  B  for  an  example  of  this  method.).  The  resulting 
stream  function  is 

=  7  )  sin  (y)  \peAx  +  qeBx  -  i),  (17) 

where  p  =  (j^L)  and  q  =  1  -  p 

In  order  to  analyze  the  stream  function  numerically,  Stommel  introduced  the  following 
parameters:  X  =  109cm,  b  =  2tt  x  108cm,  D  =  2  x  104cm,  F  =  1  dyne/cm2,  and 
R  —  0.02.  F  and  R  were  picked  arbitrarily  in  order  to  produce  the  correct  physical 
features. 

When  Coriolis  is  considered  to  be  a  linear  function  of  latitude,  the  streamlines 
indicate  an  intensification  of  current  velocities  along  the  western  boundary.  This  is  the 
region  of  the  Gulf  Stream  [Stommel,  1948]. 


2.4  Csanady’s  Transport  Stream  Function 

In  Csanady  and  Hamilton  (1988),  a  transport  stream  function  for  slopewater  circu¬ 
lation  was  introduced.  It  is  given  by 


dipdf  dipdf 
dx  dy  dy  dx 


(18) 


where  W  is  the  value  for  the  wind  stress  curl  over  the  rectangular  region.  Here,  ip  is 
the  stream  function.  The  approximations  to  the  Coriolis  acceleration  are  and 
The  transport  stream  function  is  subject  to  the  following  boundary  conditions: 


ip=  { 


0 

0 

-K^ 

(G  +  *0f 


if  x=0,  0<y<b/2 
if  y=0,  0<x<a 
if  x=0,  b/2<y<b 
if  x=a,  0<y<b 
if  y=b,  0<x<a 


(19) 


where  G  is  the  inflow  from  the  Gulf  Stream,  K  is  the  inflow  of  the  Coastal  Labrador 
Sea  Water,  x  and  y  are  the  coordinates,  a  is  the  length  of  the  box  along  the  x-axis,  and 
b  is  the  length  of  the  box  along  the  y-axis. 

The  differential  equation  of  Csanady  and  Hamilton  was  developed  using  several  impor¬ 
tant  assumptions  that  are  mentioned  below: 


•  the  dimensions  of  the  idealized  basin  are  200  km  wide  by  1600  km  long 

•  the  mean  depth  is  500  m  overlying  an  inert,  deep,  very  heavy  water  mass 

•  the  Coriolis  parameter  is  /  =  10-4s_1  which  increases  northward  at  the  rate  of 
/?  =  1.6  x  10-nm_1s-1, 
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•  the  wind  stress  curl  is  one  of  the  driving  factors  of  the  circulation, 

•  the  flow  is  steady-state,  and 

•  a  rotated  coordinate  system  is  used  where  the  cartesian  coordinate  system  is 
rotated  67  deg  from  north. 

By  inspection,  one  can  find  the  differences  between  the  equation  stated  by  Csanady 
and  Hamilton  and  the  equation  given  by  Stommel  in  the  previous  section.  Csanady  and 
Hamilton  discount  any  frictional  dissipative  force  by  discarding  the  Laplacian  term. 
This  helps  to  simplify  the  problem  and  the  numerics  associated  with  it.  The  bound¬ 
ary  conditions  are  non-zero  indicating  an  influence  imposed  upon  the  body  by  cur¬ 
rent  flows.  The  wind  stress  curl  is  a  constant  independent  of  the  y  variable.  (See 
[Csanady  and  Hamilton,  1988]  for  more  details.) 


3  Discretization 


.  The  first  steP  jn  developing  the  numerical  scheme  is  to  discretize  a  differentia]  equa¬ 
tion  which  describes  the  flow.  The  equation  selected  was  a  compilation  of  Stommel’s 
streamfunction  and  Csanady  and  Hamilton’s  stream  function.  The  equation  introduced 
in  the  previous  section  was  modified  to  include  a  frictional  dissipative  force  with  the  ad¬ 
dition  of  the  Laplacian.  The  differential  equation  that  was  the  basis  for  all  calculations 

IS 

zV2!};  -f  aipx  =  W  (20) 

where  e  and  a  are  arbitrary  constants  and  W  is  the  wind  stress  curl. 


3.1  Rotation  of  Coordinates 


The  previous  equation  is  oriented  in  the  North-South  coordinate  system.  -To  make 
the  equation  applicable  to  the  slope  sea  region  and  for  better  comparison  to  Csanady 
and  Hamilton  s  equation,  it  is  necessary  to  rotate  the  coordinates  by  67  deg  from  North 
The  rotation  was  accomplished  in  the  following  linear.  In  general- 


where 

and 


/ 

r 

r 

x' 

/ 

— 

cos  a 

sin  a 

X 

¥ 

—  sin  a 

cos  a 

y 

x  cos  a  +  y  sin  a, 


yf  —  — x  sin  a  +  y  cos  a 

are  the  rotated  coordinates.  In  the  slopewater  case: 


cos(27t  -  a) 
-sin(2ir  -  a) 


sin(2x  -  a) 
cos(2 7r  -  a) 


(21) 

(22) 

(23) 


(24) 
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Using  trigonometric  identities,  the  previous  equation  becomes 


xf  cos  a  —  sin  a  x 

y1  sin  a  cos  a  y 

(25) 

Here  the  rotated  coordinates  are  given  by 

x  —  x  cos  a  —  y  sin  a, 

(26) 

y*  =  x  sin  a  +  y  cos  a 

(27) 

By  utilizing  the  previous  two  identities,  it  is  possible  to  rotate  Stommel’s  equation 
into  Csanady  and  Hamilton’s  region  for  Slopewater.  By  performing  this  operation,  it 
becomes  apparent  how  the  slopewater  transport  function  was  developed.  The  mathe¬ 
matical  operations  are  outlined  below: 

Wx,y)  =  ^'(x\yr) 

(28) 

,  ,  ^  dil)' dx'  dip'  dy' 

*^’y)=a?dI+dy'd, 

(29) 

dip'  dip'  . 

ipx  =  ~Q^j cosa  +  —  sin  a 

(30) 

,  d  (dip'  dip'  .  \ 

'/’™  =  te[#?cosa+ws,nV 

(31) 

Let  cos  a  =  c  and  sin  a  =  s,  then  the  second  derivatives  of  the  streamfunction  become 


,  j  d2ip'  ,  n  d2ip'  ,  2  d2ip' 

=  c  dx°  +  2cs  dx'dy' +  s  ay * 

(32) 

(  dip'  dx'  dip'  dy' 

v  dx'  dy  '  dy'  dy 

(33) 

dip'  dip’ 

^  =  ~d?s  +  Wc 

(34) 

,  2  d2ip'  d2ip'  2d2ip' 

*”  *  c  B, «  2csdx'dy'  +  5  dx* 

(35) 

The  Laplacian  V2^  is  defined  to  be  +  ipyy.  When  this  is  performed  on  the  equa¬ 
tions  for  rotated  coordinates,  the  sine  and  cosine  terms  drop  out  leaving  the  following 
equation: 

$ XX  +  ifcyy  =  Vx'x'^y'y'  (36) 

From  this  point  onward,  all  references  to  x,y,and  'ip  will  be  assumed  to  be  in  the  rotated 
coordinate  system  unless  stated  otherwise. 

The  new  equation  in  the  rotated  coordinate  system  becomes  by  substitution  of  the 
previous  equation  into  Csanady  and  Hamilton’s  transport  streamfunction 

^(^xx  +  tyyy)  +  Cx(c7px  +  S'tpy)  =  W 


(37) 
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3.2  Finite  Differences 

Finite  differencing  is  a  numerical  method  used  to  describe  a  continuous  region  with 
discrete  points.  At  each  point  there  is  an  approximated  value  which  describes  the  do¬ 
main.  A  truncation  error  is  produced  for  each  approximation.  This  will  be  the  difference 
between  the  exact  solution  and  the  numerical  scheme.  The  method  of  approximation 
is  based  on  the  expansion  of  the  Taylor  series.  For  centered  differences,  the  Taylor 
senes  is  expanded  both  in  the  forward  direction  and  in  the  backward  direction,  i.e. 

u(X  +  Ax,  i,)  and  u(x  -  A x,y),  subtracted,  and  divided  by  twice  the  incremental  step 
(See  [Ames,  1992]  for  specific  details.)  . 

For  this  problem  the  Laplacian  of  ip  and  the  first  order  derivatives  of  ip  must  be 
approximated  using  the  following  finite  centered  differences  (Appendix  A): 

ib  ~  ~  +  tidzl  /OON 

%y  -  p  (38) 


i’xx  ^ 
tPy  ~ 
lbx  ~ 


y,-+ hi  -  flid  ±  Yi-Li 


h 2 

2  k 

ti+hi  ~  tizld 
2  h 


(39) 

(40) 

(41) 


where  h  _  ,  k  -  and  i  and  j  reference  the  discrete  points  within  the  domain. 

These  appr^mations  are  substituted  into  the  transport  streamfunction  equation 
with  rotated  coordinates.  The  finite  differenced  equation  becomes 

£  +  ,  j  + 

O,  [c  +  ^  j  =  ^ 

Solving  for  one  is  left  with  the  approximation  of  the  domain  at  each  point. 

ib-  ■  =  —  -  .  a3  /  a4  a5 

t,:  al  al  ‘,J+1  al^,,J-1  ~  al^,+1’'7  ~  (42) 

where  al  = -2£  (^  +  J,) ,  a2  =  S,  „3  =  ^ and  „5  =  ^  -  g. 

3.3  Iterative  Method 

The  iterative  method  used  for  the  numerical  scheme  is  the  Gauss-Seidel  method  This 
method  owes  its  derivation  from  the  linear  algebra  equation  that  is  given  by  Ax  =  b. 
This  method  is  known  as  the  method  of  successive  displacements  or  iteration  by  single 
s  eps  ere,  the  calculations  are  based  on  the  immediate  use  of  the  improved  values 
(  ee  [Ames,  1992]  for  more  details.).  The  computation  utilizes  the  following  equation: 
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(fc)  _  J2)=\(aijx\  -  EjU'+l  (a*'jxj  1)  +  ^t) 

**’  ”  "  —  (43) 

uu 

where  k  is  the  number  of  iterations  and  i  and  j  correspond  to  the  row  and  column  entries 
in  the  matrix. 

In  the  numerical  scheme  for  slopewater  circulation  the  number  of  iterations  were 
determined  by  the  numerical  output.  The  program,  written  using  the  Mathematica 
package  and  run  on  a  SPARC20  Sun  station,  was  instructed  to  continue  computation 
until  the  maximum  error  was  on  the  order  of  10-7. 

4  Calibration  Techniques 

Numerical  schemes  must  be  validated  before  one  can  feel  confident  that  the  output  of 
the  model  is  indeed  reasonable.  This  was  accomplished  by  two  methods  of  calibration. 
The  first  was  a  comparison  to  the  Stommel  Gulf  Stream  model.  The  second  was  a 
boundary  condition  verification. 

4.1  Slopewater  vs.  Gulf  Stream 

The  first  calibration  was  a  comparison  of  numerical  output  to  an  exact  solution. (See 
Figure  2.)  The  exact  solution  to  Stommel ,s  differential  equation  discussed  previously 
was  produced.  This  was  then  compared  to  the  numerical  output  when  the  parameters  of 
Stommel ’s  model  were  introduced  into  the  model.  In  addition  the  boundary  conditions 
were  set  equal  to  zero  and  the  forcing  term  was  set  equal  to  a  sine  function  instead  of 
a  constant.  The  absolute  value  of  the  difference  between  the  maximum  values  of  the 
exact  solution  and  the  numerical  solution  was  then  calculated.  The  difference  was  on 
the  order  of  10-8  after  500  iterations. 

4.2  Boundary  Verification 

The  second  calibration  was  a  boundary  condition  verification.  The  stream  function 
was  determined  to  be  ip(x,  y )  =  sin(^)  sinh(^)  by  separation  of  variables.  (See  Appen¬ 
dices  B  and  C.)  The  first  normal  mode  was  given  as  the  northern  boundary  condition 
with  all  other  boundaries  being  zero.  In  the  numerical  scheme,  W,  a,  and  s  were  set 
equal  to  zero  while  c  was  set  equal  to  one.  The  numerical  solution  was  once  again 
compared  to  the  exact  solution  (Figure  3). 

5  Discussion  of  Model  Outputs 

Several  runs  of  the  model  were  performed  using  various  parameters.  The  first  set 
of  parameters  was  based  solely  on  the  constants  and  boundary  conditions  given  by 
Csanady  and  Hamilton.  The  second  set  of  parameters  was  a  combination  of  Csanady 
and  Hamilton’s  parameters  and  self-imposed  boundary  conditions. 


£feUL2;  °*> 

the  streamfunction.  The  axes  represent  the  dimensions  of  the  calTullfm^J^  * 


Figure  3:  Boundary  Verification  by  comparison  of  the  exact  solution  of  a  specified 
streamfunction  to  the  numerical  solution  of  the  same  streamfunction.  The  curves  are 
streamlines  with  different  shades  indicating  different  values  of  i/;.  The  exact  contour  has 
axes  represented  in  units  of  length.  The  numerical  solution  has  axes  which  represent 
the  dimensions  of  the  matrix.(See  Appendix  D.l  .) 
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5.1  Csanady  and  Hamilton’s  Parameters 

Three  different  products  were  produced  by  varying  the  Gulf  Stream  constant,  G, 
which  is  a  value  of  flux  rate  of  the  Gulf  Stream.  In  the  case  where  G  =  2  x  106m3s-1, 
there  is  an  intense  crowding  of  stream  lines  near  the  northwestern  boundary  indicating 
a  large  velocity  gradient.  There  is  a  small  deflection  eastward  of  the  Gulf  Stream.  (  See 
Figure  4.)  As  G  increases  to  6  x  106  and  lOVr1  there  is  an  increased  deflection  of 
the  Gulf  Stream  eastward  and  a  decreasing  velocity  gradient  in  the  northwest  corner. 
(See  Figures  5  and  6.).  These  products  can  be  compared  to  the  outputs  of  Csanady 
and  Hamilton.  A  comparison  shows  that  both  models  indicate  greater  deflection  of  the 
Gulf  Stream  with  increased  flux  rate. 

The  boundary  conditions  imposed  do  not  create  a  velocity  field  which  corresponds 
to  the  field  that  is  actually  observed.  (  See  Figure  1.)  The  velocity  field  only  represents 
a  small  section  of  the  slopewater  region.  In  order  to  more  accurately  represent  the 
circulation  which  is  known  to  occur,  new  boundary  conditions  were  developed  and 
introduced  into  the  model. 

5.2  New  Boundary  Conditions 

New  boundary  conditions  were  developed  from  inspection  of  the  observational  scheme. 
Based  on  the  actions  of  the  streamlines,  functions  were  chosen  to  approximate  the 
motion  of  the  flow.  The  boundary  condition  on  the  continental  shelf  was  determined 
to  be J>(0,y)  =  0  because  there  is  no  fluid  flow  across  the  boundary,  only  parallel  to 
the  shelf.  The  boundary  condition  on  the  southwestern  edge  ^(x,0)  =  0  because  the 
flow  changes  directions  but  once  again  does  not  cross  the  boundary  of  the  region.  The 
final  two  boundaries  were  approximated  by  functions  dependent  on  either  x  or  y.  The 
Gulf  Stream  side  has  two  boundary  conditions.  For  ip(l,y),  when  0  <  y  <  6/A,  ip  =  0. 
When  6/A  <  y  <  6/2A,  ip  is  approximated  by  a  line.  The  Gulf  Stream  is  represented 
in  this  fashion  because  after  the  current  flows  along  the  region  for  half  of  the  distance 
the  stream  enters  into  the  slopewater.  The  northern  edge  must  represent  the  change 
in  direction  of  the  fluid  flow  due  to  the  exit  of  the  Gulf  Stream  and  the  entrance  of 

the  Coastal  Labrador  Sea  Water.  This  physical  process  is  represented  by  a  quadratic 
equation.  (See  Appendix  C.) 

The  first  run  of  the  new  boundary  conditions  (See  Figure  7.)  showed  the  development 
of  a  western  gyre  and  the  entrance  of  the  Gulf  Stream  into  the  region.  Because  the 
output  did  not  closely  resemble  the  schematic  (See  Figure  1.),  constants  then  had  to  be 

varied  until  the  combination  was  discovered  which  would  more  accurately  represent  the 
schematic  (Figure  7). 


6  Data  Analysis 

To  test  the  credibility  of  the  numerical  calculations,  model  outputs  must  be  compared 
to  real  numbers.  Data  for  comparison  comes  directly  from  the  MASAR  experiments 
Charts  of  the  region  showing  current  meter  deployment  contain  velocity  vectors  for 
the  upper  level  currents.  To  retrieve  the  data,  the  resultant  speed  was  determined  by 
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Figure  4:  Slopewater  model  output  (See  Appendix  D.4.)  when  using  Csanady  and 
Hamilton’s  parameters  and  boundary  conditions.  In  this  case,  G  is  2  x  10^m?s  .  The 

curves  are  streamlines.  The  orientation  of  the  output  is:  the  northern  edge  is  to  the 
right  of  the  page  and  the  Gulf  Stream  is  to  the  bottom.  This  orientation  is  the  same 
for  the  next  two  products.  The  axes  indicate  the  dimensions  of  the  calculated  matrix. 
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Figure  5:  Slopewater  output  when  using  Csanady  and  Hamilton’s  parameters  and 
boundary  conditions.  In  this  case,  G  is  6  x  106m3s-1  with  the  curves  indicating  the 
streamlines.  The  axes  indicate  the  dimensions  of  the  calculated  matrix. 
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Figure  6:  Slopewater  output  when  using  Csanady  and  Hamilton’s  parameters  and 
boundary  conditions.  In  this  case,  G  is  107ro3s-1.  The  level  curves  are  streamlines. 
The  axes  indicate  the  dimensions  of  the  calculated  matrix. 
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Figure  7:  Slopewarsr  ~~ — — - ,  , 

curves  are  streamline  *Ith  new  boundary  conditions  imposed.  The  level 

orientation  of  the  ah°kW  the  number  of  nodes  used  ^  the  matrix.  The 

Stream  is  to  the  brram  --  -T'  edSe  15  to  tbe  right  of  the  page  and  the  Gulf 

closed  gyre  in  the  ^  Some  important  features  in  the  output  are  the 

Stream  into  the  re-cnT^^  ^  ?*  SL°pe  regi°n  and  the  intrusion  of  the  Gulf 

similar  but  the  Cc«--  -  *ompared  to  fi§ure  1.  one  notices  that  the  features  are 

Appendix  D  4)  Current  is  overpowered  by  the  Gulf  Stream.  (See 
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Figure  8:  Slopewater  circulation  using  new  boundary  conditions.  This  output  is  the 
result  of  many  runs  looking  for  the  coefficients  that  would  produce  features  that  closely 
resemble  Figure  1.  Note  the  closed  gyre,  the  intrusion  of  the  Gulf  Stream,  the  Coastal 
Labrador  Current,  and  what  appears  to  be  the  change  in  direction  of  the  Coastal 
Labrador  Current.  The  axes  indicate  the  number  of  nodes.  The  level  curves  are  stream¬ 
lines.  The  orientation  of  the  box  is  as  follows:  the  northern  edge  is  to  the  right  of  the 
page  and  the  Gulf  Stream  is  to  the  bottom  of  the  page. 
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measuring  velocities  on  a  chart  within  the  MASAR  experiment  report  using  the  given 
scale.  In  order  to  compare  the  velocity  vectors  to  the  stream  function,  the  resultant 
vectors  had  to  be  broken  down  into  the  u  (the  velocity  in  the  x  direction)  and  v  (the 
velocity  in  the  y  direction)  components.  The  stream  function  must  be  differentiated  to 
derive  the  numerically  calculated  u  and  v  components. 


7  Naval  Applications 

The  understanding  of  the  coastal  environment  is  important  for  military  operations 
in  the  littoral  environment.  The  dynamics  within  the  area  of  operation  affect  many 
different  warfares.  In  particular,  knowledge  of  currents  can  be  beneficial  for  successful 
mine  operations,  submarine  warfare  and  anti-submarine  warfare.  Currents  can  even  aid 
or  hinder  navigation. 

The  slopewater  region  is  potentially  a  very  important  tactical  region.  The  entire 
eastern  seaboard  is  bounded  by  the  water  mass.  The  Gulf  Stream  has  a  large  influence  in 
the  region  and  warm  core  eddies  are  typically  spun  off  in  this  area.  These  eddies  and  the 
fast  moving  current  of  the  Gulf  Stream,  if  not  fully  understood,  could  hinder  navigation 
and  mine  sweeping.  The  change  in  the  temperature  profile  alters  the  sound  velocity 
profiles  when  warm  core  eddies  are  interacting  with  the  slopewater  mass.  This  provides 
submarines  a  better  chance  to  remain  undetected.  Knowledge  and  understanding  of  the 
circulation  of  the  slopewater  and  how  the  eddies  move  within  the  circulation  should  give 
the  tactical  commander  a  better  idea  of  where  to  search  for  the  enemy  or  where  to  hide. 
The  currents  in  the  slopewater  are  not  as  fast  along  the  shelf  as  they  are  along  the  Gulf 
stream  region.  Understanding  this  and  being  able  to  determine  the  speed  changes  as 
you  move  shoreward  will  benefit  navigation. 

The  modeling  of  slopewater  is  applicable  to  many  gyres  throughout  the  world’s 
oceans.  The  same  techniques  can  be  used  to  generate  numerical  schemes  for  any  type 
of  wind-driven  circulation.  If  the  Navy  is  interested  in  better  understanding  the  coastal 
environments,  numerical  modeling  is  another  approach  that  can  be  utilized. 


8  Future  Developments 

The  numerical  scheme  developed  is  highly  simplified.  It  is  two-dimensional,  time 
independent,  and  based  on  a  limited  number  of  physical  parameters.  The  region  is 
rectangular  and  considered  to  be  homogeneous  and  to  have  a  constant  depth.  The  wind 
stress  curl  is  constant  and  the  flux  rate  of  the  Gulf  Stream  is  constant.  The  dynamic 

implications  of  warm  core  eddies  and  changes  in  temperature,  pressure,  and  salinity  are 
ignored. 

In  order  to  improve  the  model,  each  of  these  factors  should  be  addressed.  It  is  not 
a  reasonable  expectation  to  have  every  influence  play  an  accurate  role  in  the  modeling 
of  the  region.  At  the  present  time,  the  knowledge  and  ability  is  unavailable  to  tackle 
a  problem  of  that  magnitude.  The  parameters  must  be  addressed  either  one  or  two  at 
a  time.  For  example,  the  region  could  have  non-rectangular  dimensions  developed  and 
the  wind  stress  curl  could  be  made  into  a  function  instead  of  a  constant.  There  are 
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many  avenues  of  approach  depending  on  the  parameters  of  interest. 

Other  improvements  could  be  made  solely  on  the  mathematical  techniques  of  com¬ 
putation.  Speed  of  computation  and  truncation  error  can  be  improved  in  a  number  of 
ways.  The  discretization  method  could  be  altered  to  produce  a  non-uniform  grid  which 
allows  for  better  mapping  of  the  region  of  interest.  The  iterative  technique  could  be 
changed  as  well  in  order  to  increase  speed  of  computation.  There  are  several  avenues  of 
approach  to  numerical  improvements  which  develop  in  higher  level  mathematics. 

The  numerical  scheme  has  considerable  room  for  improvement.  However,  the  work 
accomplished  is  not  trivial.  The  product  resembles  what  is  generally  understood  to  be 
true  about  the  circulation  in  the  region.  Improvements  would  allow  for  more  accuracy 
and  variation  in  the  model  outputs. 


24 


References 

Ames, William  F  (1992).  Numerical  Methods  for  Partial  Differential  Equations.  Aca¬ 
demic  Press,  Inc.,  New  York,  451pp. 

Csanady,  G.T.  and  P.  Hamilton  (1988).  Circulation  of  Slopewater.  Continental  Shelf 
Research ,  8,  565-624. 

Mid-Atlantic  Continental  Slope  and  Rise  Final  Report-  May  1987,  2. 

Munk,  W.H.  (1950).  On  the  wind-driven  ocean  circulation.  Journal  of  Meteorology,  7, 
79-93. 


Neumann,  G.  (1968).  Ocean  Currents.  Elsevier  Scientific  Publishing  Company,  New 
York,  352pp. 

0  Neill,  M.E.  and  F.  Chorlton  (1986).  Ideal  and  Incompressible  Fluid  Dynamics.  Ellis 
Horwood  Limited,  New  York,  412pp. 

Pedlosky,  J.  (1987).  Geophysical  Fluid  Dynamics.  Springer- Verlag,  New  York,  710pp. 

Pond,  S.  and  G.L.  Pickard  (1983).  Introductory  Dynamical  Oceanography.  Pergamon 
Press,  New  York,  329pp. 

Stommel,  H.  (1958).  The  Gulf  Stream.  University  of  California  Press,  Los  Angeles, 

202pp. 


Stommel,  H.  (1948).  The  westward  intensification  of  wind-driven  ocean  currents.  Trans¬ 
actions,  American  Geophysical  Union ,  29,  202-206. 


A  Development  of  Finite  Difference  Equations 
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Add  the  four  equations  together  to  get  the  following: 


V’.+ij  +  V’i-i.j  +  V’i.j+i  +  V’t.j-i  —  4'0t,j  +  h2tpxx  +  k  tpyy 

Therefore, 

V’xx  ~  -  2 V’.J  +  V’.-I.i) 

’’/’yy  ~  £2  (V’t.j+l  —  2 ipij  +  V’t.j-l) 

Given 


V’x  ~  -  V’t-i.j) 

V>y  ~  -  V’.J-l) 


/l  = 


A 

n  +  1 


m  +  1 

Substitute  into  the  PDE  for  slopewater  circulation. 


e{lpxx  +  t/’yy)  +  a(cV>x  +  StM  =  W 
£  +  iK+i., -2^,,  +0-1,^  + 


,+  1  ~ V'i.j-l 

2k 


ahpij  =  W  -  a2ipitj+\  -  a3V>;j-i  -  a40,+i,j  -  o5^«-ij 
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where  al  =  —2s  (p  +  p) 

«2=^r+ff 
a3  =  pr  “  ff 
«4  =  ^+f 

a5  =  ^r  +  f^ 


B  Separation  of  Variables 


^Pxx  “i”  ^Pyy  —  0 

Boundary  Conditions 

^(0 ,y)  =  0 
V>(x,0)  =  0 
V»(A  ,y)  =  0 
ij)(x,b)  =  0 

Suppose 

=  F{x)G{y) 

Then 

i’xx  =  F"G 
i>yy  =  FG" 

Substitute 

F"G  =  -FG" 

Divide  by  F"G" 

G_  _  _ 

G"  ~  F"  ~ 

G  =  -G"fi 
G  +  G"n  =  0 

let  y/^JI  =  k 


G{y)~  ci  sinh  ky  +  c2  cosh  ky 
-F  =  -nF" 

HF"  -  F  =  0 


F{x)  =  c3  sin  kx  +  c4  cos  kx 
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The  general  solution  of  y>  becomes 

rp(x,  y )  =  (ci  sinh  ky  +  c2  cosh  ky)(c3  sin  kx  +  c4  cos  kx ) 
To  solve  for  the  exact  solution,  boundary  conditions  must  be  applied. 

V»(0  ,y)  =  0 


0  =  (ci  sinh  ky  +  c2  cosh  ky)(c3  sin  0  +  c4  cos  0) 


0  =  C4(ci  sinh  ky  +  c2  cosh  k ) 


^>(x,0)  =  0 


c4  =  0 


0  =  (ci  sinh  0  +  c2  cosh  0)(c3  sin  kx) 


0  =  (c2)(c3sin  kx) 
c2  =  0 


ip(x,  y)  =  (c!  sinh  ky)(c3 sin  kx) 

Let  C1C3  =  c  =  1 
i>(\,y)  =  0 


0  =  sinh  ky  sin  kX 


0  =  sin  k  A 


kX  =  717T 


Therefore  the  exact  equation  is 


f/  N  .  titcx  .  ,  niry 
rp{x,  y)  =  sin  — -  smh  — — 

A.  A 

The  boundary  condition  for  tp(x,b)  is  designated  to  be  the  first  normal  mode  of  the 
exact  solution.  It  is  given  by 


4>(x,  b)  =  sin  ^  s*n^  ~T 
A  A 
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C  Development  of  Boundary  Conditions 

The  Shelf  boundary  condition  is  ip(0,y)  =  0. 

The  Southern  edge  boundary  cpndition  is  ip(x,  0)  =  0. 

The  Gulf  Stream  boundary  condition  is  ip(l,y)  =  0  when  0  <  y  < 

Where  £  <  V  <  I  »  Why)  =  *p(v-  &) 

^  is  a  point  through  which  the  line  must  pass.  This  point  represents  the  maximum 
flux  of  the  Gulf  Stream. 

The  Gulf  Stream  boundary  orginates  from  the  equation  of  a  line  which  passes 
through  the  points  (jy,0)  and  (%,£).  The  derivation  is  as  follows: 

y  =  mx  +  B 


y-B 

m 


2X6  b  , 


6(£+£2)-(r?+6£2 

£+F 


x2  + 


The  Northern  edge  boundary  condition  is  i/>(x,j)  = 

This  boundary  condition  is  based  on  a  quadratic  equation  which  must  pass  through  the 
points  (0,0),  (£,-»?),  and  (1,£).  The  derivation  is  as  follows: 


f(x)  —  ax2  +  bx  +  c,  where  a,b,  and  c  are  constants  to  be  determined  below. 

To  solve  for  the  constants  a,b,  and  c,  substitute  the  coordinates  into  the  quadratic 
equation. 

Solve  for  c. 


at  (0,0) 


0  =  0  +  0  +  c 

at  (f ,  -77) 

c  =  0 

at  (1,6) 

-Tj  =  of2  +  b£ 

Therefore,  a  =  6  —  b. 

Solve  for  b. 

6  =  a  +  b 

-T)  =  a£2  +  added  to  -6£2  =  -a£2  -  b£2  yields 


»7  +  *£2  =  &(f  +  f2) 
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Solve  for  a  by  substitution. 


_  v±se 
t+e 


So, 


a  = 


^+e2)-o?+^2) 


1p(x,j)  = 


+  t2)  -  (V  +  6?) 


t+e 


2  , 

rj  +  6£2 

X  + 

e  +  e. 

x 


D  Mathematica  Programs 

D.l  Boundary  Verification 

(*  dimensions  of  box  *) 

lambda=10"9; 

b=2*Pi*10"8; 

(♦depth  of  layer*) 
d=2*10‘4; 

(♦force  of  wind*) 

F  =  1; 

(♦frictional  dissipative  term*) 

R=0.02; 

(*  constants  *) 
gamma=F*Pi/R/b; 
eta=N[Pi/b] ; 
coriolis=10~(-13) ; 

(*  dimensions  of  matrix  *) 
m=  64; 
n=  64; 

(*  incremental  step  size  *) 
k=b/lambda/(m+l) ; 
h=l/(n+l) ; 

disx=Table[N[i*h]  ,-Ci,l,n>]  ; 


disy=Table [N [j  *k] , { j , 1 ,m>] ; 


(*  constants  *) 


eps=  1; 
alpha=  0; 

(*  cosine  *) 
c=0; 

(*  sine  *) 
s=0 ; 

(*  wind  stress  curl  *) 

W=Table[0,{i,m+2>]  ; 

(*  constants  *) 

al=  N[-2*eps*(l/(k-2)+l/(h‘2))]; 

a2=  N[eps/(k~2)+alpha*s/(2*k)] ; 

a3=  N[eps/(k‘2)-alpha*s/(2*k)] ; 

a4=  N[eps/(h~2)+  alpha*c/(2*h)] ; 

a5=  N[eps/(h~2)-  alpha*c/ (2*h)] ; 

al2  =  a2/al;  al3=a3/al;  al4=a4/al;  al5=a5/al; 

(*  Boundary  Conditions  *) 

gl[y_]  =  0;  fl[x_]  =  0;  g2[y_]  =  0; 
f2[x_]=  N  [Sin  [Pi  x]  *Sinh[Pi]]  ; 

psi=Table[0,{i,m+2},{j ,n+2>]  ; 
newpsi  =  psi; 

Do[psi[[l,  j  +  l]]=gl[disy[[j]]]  ,  { j  ,  m>]  ; 
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Do[psi[[n+2, j+l]]=  g2[disy[[j]]]  ,  {j ,  m>]  ; 

Do[psi[[i+l,l]]=  f 1 [disx[[i]]] ,  {i,  n>] ; 

Do[psi[[i+l,m+2]]=  f2[disx[[i]]]  ,  {i,  n>]  ; 

Do  [ 

Do [newpsi [ [i , j] ] =  W[[i-l]]/al-al2*psi[[i, j+1]] -al3*psi [[i , j-l]] 
-al4*psi[[i+l , j]] -al5*psi [[i-1 , j]] ,  {i,2,n+l>,  {j ,2,m+l>] J 
error=Max [Abs [psi-newpsi]] ;  Print [error] ; 

Do[psi [[i,  j]]=newpsi[[i, j]] ,  {i,  2,  n+l>,  { j  ,  2,  m+l>] , 

-Cii,  50>] 


AA=ListContourPlot [Transpose [psi] ,ColorFunction->Hue] 

exactsoln [x_ ,y_]=Sin [Pi  x]*Sinh[Pi  y*lambda/b] 

zz=Table [exactsoln [disx [ [i] ] ,disy[[j]]]  ,{i,n+2>,{j ,m+2>] 

ZZ=ContourPlot  [exactsoln [x , y] , {x , 0 , 1} , {y , 0 , b/lambda} , 
ColorFunction->  Hue] 

AAA=Show [AA ,PlotLabel->Mnumerical"] 

PSPrint [AAA] 

zzz=Shov[ZZ,PlotLabel->"  exact"] 

PSPrint [zzz] 

Max [Abs [zz-newps i] ] 

AAAA=  Show  [Graphics Array  [{{AAA} ,  {zzz}}]  ] 

Display ["graphl.ps",  AAAA] 


PSPrint [AAAA] 


D.2  Stommel  Verification 


(*  dimensions  of  box  *) 
lambda=2*10~7 ; 
b=16*10"7 ; 

(*  depth  of  layer  *) 
d=2*10~4; 

(*  force  of  wind  *) 

F  =  1; 

(*  frictional  dissipative  term  *) 

R=0.02; 

(*  constants  *) 
gamma=F*Pi/R/b*lambda ; 
eta=N[Pi/b]  ; 
coriolis=10“(-13) ; 

(*  Gulf  Stream  flux  rate  *) 

G=2*10‘12; 

(*  Coastal  Labrador  Sea  Water  flux  rate  *) 
K=4*10~12; 

(*  matrix  dimensions  *) 
m=  64; 
n=  64; 


(*  incremental  step  size  *) 
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k=b/lambda/ (m+1) ; 
h=l/ (n+1) ; 

disy=Table [N [i*k] ,{i ,  1 ,m>] ; 
disx=Table [N  [- j  *h] , { j , 1 , n}] ; 

(*  constants  *) 
eps=  1; 

alpha®  l*lambda; 

(*  cosine  *) 
c=0; 

(*  sine  *) 
s=  0; 

(*  wind  stress  curl  *) 

W=Table[l,<i,l,m>]; 

C*  constants  *) 

al=  N[-2*eps*(l/(k-2)+l/(h‘2))]  ; 

a2=  N[eps/(k“2)+alpha*s/(2*k)] ; 

a3=  N[eps/(k‘2)-alpha*s/(2*k)]  ; 

a4=  N[eps/(h~2)+  alpha*c/ (2*h)] ; 

a5=  N[eps/(h“2)~  alpha*c/(2*h)]  ; 

al2  =  a2/al;  al3=a3/al;  al4=a4/al;  al5=a5/al; 

(*  Boundary  Conditions  *) 

gl [x_] =0 ; g2 [x J  =0 ; f 1 [y J  =0 ; f 2 [y J  =0 ; 
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psi=Table[0,{i,m+2},-Cj,n+2}]  ; 
newpsi  =  psi; 

Do  [psi  [[1 ,  j+l]]=gi  [disx [[j]]]  ,  -Cj,  n}]  ; 

Do[psi[[m+2,  j+l]]  =  g2[disx[[j]]]  ,  { j  ,  n>]  ; 

Do[psi[[i+l,l]]=  f l[disy[[i]]] ,  {i,  m>] ; 

Do [psi [ [i+1 ,n+2]] =  f2[disy[[i]]] ,  {i,  m>] ; 

Do[ 

Do [newpsi H[[i-l]]/al-al2*psi[[i, j+l]]-al3*psi[[i, j-1]] 
-al4*psi[[i+l,j]]-al5*psi[[i-l,j]]  ,  -[i,2,m+l},  -Cjf2,n+1}]; 
error=Max [Abs  [psi-newpsi]] ;  Print  [error] ; 

Do[psi[[i,  j]]=newpsi[[i,j]] ,  {i,  2,  m+l>,  { j ,  2,  n+l>] , 

{ii,  50}] 

AA=ListContourPlot [psi,  PlotLabel->Stommel  Verification, 
ColorFunction->Hue,Contours->20] 

PSPrint [AA] 

Display["stommelver.psM ,  AA] 


D.3  Csanady’s  Boundaries 


(*  Gauss-  Seidel,  Stommel  *) 

(*  The  real  psi  is  related  to  this  by  psi  by 
lambda*  this  psi(x/lambda,  y/lambda)*) 

(*  dimensions  of  box  in  meters  *) 

lambda=2*10~5; 

b- 16*  1CT5 ;  scale=lambda/b ; 

(*  depth  *) 

d=5*lCT2; 

(*  force  of  wind  *) 

F=2. 4*1(T-10; 

(*  frictional  dissipative  term  *) 

R=0 . 0002 ; 

(*  constants  *) 

gamma=F/R ;  gammabar=gamma*lambda ; 

(*  beta  plane  approximation  for  Coriolis  *) 
beta=  1.6*10^-11; 

(*  constants  *) 

alpha=  d*beta/R;  alphabar=alpha*lambda; 

(*  incremental  step  size  *) 
k=b/lambda/(m+l) ; 
h=l/(n+l) ; 


(*  size  of  matrix  *) 
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m=  41; 
n=  41; 

disx=Table[N[i*h] ,{i,l,n>] ; 
disy=Table[N[j*k] ,l,m>] ; 

(*  constants  *) 
eps=  1; 

c=  N [Cos [1.168]] ; 
s=  N [Sin [1.168]] ; 

(*  Gulf  Stream  flux  rate  *) 

G=  2*10*6  (*m“3/s*) ; 

(*  Coastal  Labrador  Current  flux  rate  *) 

K=  4*10*6  (*m*3/s*) ; 

(*  wind  stress  curl  *) 

W=Table[gammabar,-Cj  ,l,m>]  (*m/s*) ; 

(*  constants  *) 

al=  N[-2*eps*(l/(k*2)+l/(h“2))] ; 

a2=  N[eps/(k“2)+alphabcir*s/(2*k)]  ; 

a3=  N[eps/(k‘2)-alphabcir*s/(2*k)]  ; 

a4=  N[eps/(h“2)+  alphabELT*c/(2*h)]  ; 

a5=  N[eps/(h*2)-  alphabar*c/(2*h)] ; 

al2  =  a2/al;  al3=a3/al;  al4=a4/al;  al5=a5/al; 


(*  Boxmdary  Conditions  *) 


g2[y_]  =  G*y/b; 

gl[y_]  =  lf[0  <y <b/2/lambda , 0 , -K* (2*y*lambda-b) /b/lambda]  ; 
fl[x_]  =  0; 

f2[x_]  =  (G+K)*x/lambda; 
ps i=Table [0 . 1 1 ,{i ,n+2> , {  j ,m+2>] ; 

Do [psi [ [1 , j+1] ] =gl [disy [ [j] ] ] ,  { j , • m>] ; 

Do[psi[[n+2, j+l]]=  g2[disy[[j]]] ,  { j ,  m>] ; 

Do[psi[[i+l,l]]=  f 1  [disx[[i]]] ,  {i,  n>] ; 

Do[psi[[i+l,m+2]]=  f 2 [disx [[i]]] ,  {i,  n}]  ; 
oldpsi  =  psi;  error  =1;  ii=l; 
newdisx=Append [Prepend [disx,  0],  N[l]]; 
newdisy=Append [Prepend [disy ,  0],  N[l/scale]]; 

Print ["m  =  " ,  m,  "  n  -  ",  n] ; 

While [error  >  10“-7, 

Do[psi[[i, j]]=W[[j-l]]/al-al2*psi[[i, j+l]]-al3*psi[[i, j-1]] 
-al4*psi[[i+l,j]]-al5*psi[[i-l,j]] ,  {i,2,n+l},  {j,2,m+l}]; 
error=Max[Abs[psi-oldpsi]] ;Print[ii,"  ",  error]; 
Do[oldpsi[[i,  j]]=psi[[i, j]] ,  {i,  2,  n+l>,  { j ,  2,  m+l>] ; 
ii=ii+l] ; 


AAl=ListContourPlot  [Transpose [psi] ,  PlotLabel->"psi"  ,  Contours->20, 
ColorFunction->  Hue,  AspectRatio->3,Color0utput->CMYKColor] 

PSPrint [AA1] 


Display ["graph3. ps" ,  AA1] 
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D.4  New  Boundaries 

(*  Gauss-  Seidel,  Stommel  *) 

(*  Th®  real  psi  is  related  to  this  by  psi  by 
lambda*  this  psi(x/lambda,  y/lambda)*) 

(*  dimensions  of  box  *) 

lambda=l; 

b=10  ;scale=lambda/b; 

(*  depth  *) 
d=l; 

(*  force  of  wind  *) 

F=1 ; 

(*  frictional  dissipative  term  *) 

R=l; 

(*  constants  *) 

gamma=F /R ; gammabar=gamma*lambda ; 
beta=  1; 

alpha=  d*beta/R;  alphabar=alpha*lambda; 

(*  incremental  step  size  *) 
k=b/lambda/(m+l) ; 
h=l/(n+l) ; 

(*  size  of  matrix  *) 
m=  64; 


n=  32; 


disx=Table[N[i*h] ,{i,l,n>] ; 
disy=Table [N  [j  *k] ,{ j , 1 ,m>] ; 

(*  constants  *) 
eps=  1; 

c=  N [Cos [1.168]] ; 
s=  N[Sin[l .  168]]  ; 

(*  Gulf  Stream  flux  rate  *) 

G=  1  (*m“3/s*) ; 

(*  Coastal  Labrador  current  flux  rate  *) 

K=  1  (*m“3/s*) ; 

(*  constants  *) 
mu=0 .25; 
eta=l; 

(*  vind  stress  curl  *) 

W=Table [gammabar , 1 ,m}]  (*m/s*) ; 

(*  constants  *) 

al=  N[-2*eps*(l/(k“2)+l/(h“2))]; 

a2=  N[eps/(k“2)+alphabar*s/(2*k)]  ; 

a3=  N[eps/(k“2)-alphabar*s/(2*k)]  ; 

a4=  N[eps/(h*2)+  alphabar*c/(2*h)]  ; 

a5=  N[eps/(h“2)-  alphabar*c/(2*h)] ; 

al2  =  a2/al;  al3=a3/al;  al4=a4/al;  al5=a5/al; 


(*  Boundary  Conditions  *) 

gulf [y_]  »  If [0<y<b/2/lambda,0,2*G*lambda/b*(y-b/2/lambda)] ; 
shelf  [y_]  =  0; 
south [x_]  =  0; 

north  [x_]  =  ((G*(mu+mu~2)-(eta+G*(mu“2)))/(mu+mu''2))*x"2 
+ (eta+G* (mu‘2) ) / (mu+mu*2) *x ; 

psi=Table[0.11,{i,n+2},-[j  ,m+2}]  ; 

Do[psi[[l, j+l]]=shelf [disy[[j]]] ,  -Cj  ,  m>]  ; 

Do [ps i [ [n+2 ,  j  + 1]  ]  =  gulf [disy[[j]]] ,  { j  ,  m>] ; 

Do[psi[[i+l , l]]=south[disx[[i]]] ,  {i,  n>] ; 

Do[psi[[i+l,m+2]]=  north [disx[[i]]] ,  {i,  n}] ; 

oldpsi  =  psi;  error  =1;  ii=l; 

newdisx=Append [Prepend [disx,  0],  N[l]]; 

newdisy=Append [Prepend [disy,  0],  N[l/scale]]; 

Print ["m  =  ",  m,  "  n  =  ",  n] ; 

While [error  >  10“-7, 

Do [psi [[i, j]]=W[[j-l]]/al-al2*psi[[i, j+l]]-al3*psi[[i, j-1]] - 
al4*psi[[i+l,j]]-al5*psi[[i-l,j]] ,  {i,2,n+l>,  {j,2,m+l}]; 
error=Max[Abs [psi-oldpsi]] ; Print [ii,"  ",  error]; 

Do [oldpsi [[i,  j]]=psi[[i, j]] ,  {i,  2,  n+l>,  { j ,  2,  m+l>] ; 
ii=ii+l] ; 


AAl=ListContourPlot [Transpose [psi] ,  PlotLabel->"psi  ", 
Contours->  20,  ColorFunction->  Hue,AspectRatio->3, 
ColorOutput->CMYKColor] 

PSPrint [AA1] 


Display ["graph6.ps" ,  AA1] 


